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Abstract 

We  have  used  a  novel  computer  with  a  novel  integra¬ 
tion  technique  to  study  the  evolution  of  the  entire  planetary 
system  for  nearly  100  million  years.  This  calculation  coniinns 
that  the  evolution  of  the  Solar  System  as  a  whole  is  chaotic, 
with  a  remarkably  short  timescale  of  exponential  divergence 
of  about  4  milUon  years.  Additional  numerical  experiments 
indicate  that  the  dynamical  evolution  of  the  Jovian  planets 
is  chaotic  apart  from  the  rest  of  the  system,  as  is  the  motion 
of  Pluto. 


Advances  in  computer  technology  have  made  it  possible  to  begin  to  di¬ 
rectly  address  the  age-old  question  of  the  nature  of  the  long-term  evolution  of 
the  Solar  System,  with  startling  results.  Sussman  and  Wisdom  [1]  presented 
numerical  evidence  that  the  motion  of  Pluto  is  chaotic,  with  a  timescale  for 
exponential  divergence  of  nearby  trajectories  of  only  about  20  million  years. 
Subsequently,  Laskar  [2]  found  numerical  evidence  of  the  chaotic  evolution 
of  the  whole  Solar  System  excluding  Pluto,  udth  a  timescale  for  exponential 
divergence  of  only  about  5  million  years.  Laskar*s  calculation  was  feasible  be¬ 
cause  he  analytically  averaged  the  equations  of  motion  to  remove  the  rapid 
variations  with  timescales  of  order  the  orbital  period.  The  averaged  equations 
are  perturbative  and  necessarily  truncated  after  a  particular  order  in  eccentric¬ 
ity,  inclination  and  mass  ratio.  A  new  integration  of  the  whole  Solar  System 
without  these  approximations  was  clearly  required. 

Direct  integrations  of  the  whole  planetary  system  are  computationally  ex¬ 
pensive.  Notable  long-term  integrations  of  the  outer  Solar  System  include: 
the  classic  1  million  year  integration  of  Cohen,  Hubbard,  and  Oesterwinter  [3], 
the  5  million  year  integration  of  Kinoshita  and  Nakai  [4],  the  210  million 
year  integration  performed  on  the  Dijptal  Orrery  [5],  the  100  million  year 
integration  of  the  LONGSTOP  project  [6],  and  the  845  million  year  Digi¬ 
tal  Orrery  integration  of  Sussman  and  Wisdom  [1].  Studies  of  the  long-term 
evolution  of  the  whole  Solar  System  have  been  more  limited  because  the  com¬ 
putational  resources  required  are  signiiicantly  larger,  by  about  two  orders  of 
magnitude.  Integrations  of  the  whole  Solar  System  include:  the  3  million  year 
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Digital  Orrery  integration  (excluding  Mercury)  [5],  the  2  million  year  integra¬ 
tion  of  Richardson  and  Walker  [7],  and  the  recent  ±3  million  year  integration 
of  Quinn,  Tremaine,  and  Duncan  [8,  9]  (hereafter  QTD). 

We  have  developed  new  computational  techniques  and  new  computer  hard¬ 
ware  to  make  possible  a  direct  integration  of  the  whole  Solar  System  spanning 
a  significantly  longer  interval  than  previously  achieved.  Our  new  direct  inte¬ 
gration  of  the  equations  of  motion  spans  36,000,000,000  days,  or  about  98.6 
million  years.  Our  earlier  result  concerning  the  chaotic  motion  of  Pluto,  as  well 
as  the  result  of  Laskar  that  the  Solar  System  is  chaotic  are  both  confirmed.  In 
order  to  localize  the  sources  of  the  chaotic  behavior  we  have  carried  out  nu¬ 
merous  additional  long-term  integrations.  We  have  found  that  the  evolution 
of  the  Jovian  planets  is  independently  chaotic,  as  is  the  motion  of  Pluto. 

Method  of  Integration 

We  use  the  symplectic  n-body  mapping  method  of  Wisdom  and  Holman  [10] 
to  integrate  the  planetary  system.  This  new  mapping  method  is  a  generaliza¬ 
tion  of  the  mapping  method  of  Wisdom  [11,  12].  The  basic  idea  is  as  follows. 

The  Hamiltonian  for  the  planetary  n-body  problem  can  be  written 

B  —  -HjCcpter  "b  Bjntgraetienj  (1) 

where  the  first  term  represents  the  Keplerian  motion  of  each  of  the  planets  with 
respect  to  the  sun,  and  the  second  term  describes  the  planetary  perturbations. 
The  amplest  form  of  the  map  is  obtained  by  adding  short  period  terms  so  that 
the  Hamiltonian  becomes 


Bfttap  —  B^epUr  "I"  ■Hffi*eroct»«»»29r52ir(Dt),  (2) 

where  ^3«(t)  is  a  periodic  sequence  of  Dirac  delta  functions  with  period  2z’, 
and  n  is  the  mapping  frequency.  The  availing  principle  [13]  suggests  that  the 
additional  short  period  terms  do  not  significantly  aifect  the  long-period  evolu¬ 
tion.  This  Hamiltonian  is  locally  integrable:  between  the  delta  functions  the 
motion  is  purely  Keplerian  and  integration  across  the  delta  functions  is  easily 
carried  out  in  Cartesian  coordinates.  The  resulting  map  of  the  phase  space 
onto  itself  is  a  composition  of  canonical  transformations  and  is  thus  canonical 
(i.e.,  the  Jacobian  is  a  symplectic  matrix).  In  general,  to  make  a  mapping  of 
this  sort  each  part  of  the  Hamiltonian,  considered  by  itself,  must  be  not  only 
integrable  but  also  efficiently  solvable.  The  Keplerian  phase  of  the  evolution 
can  be  efficiently  carried  out  directly  in  Cartesian  coordinates  with  the  aid  of 
Gauss’s  /  and  g  functions.  The  uniform  use  of  the  Cartesian  coordinates  also 
eliminates  any  need  for  costly  intermediate  coordinate  transformations. 

The  mapping  method  can  be  made  accurate  to  arbitrarily  high-order  in  the 
mapping  step-size.  This  and  other  refinements  are  discussed  in  Wisdom  and 
Holman  [10].  In  the  applications  presented  here  the  second  order  version  of  the 
n-body  mapping  is  used.  Second  order  is  achieved  by  evolving  the  system  with 
the  Kepler  Hamiltonian  for  a  half  mapping  step,  followed  by  an  alternating 
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succession  of  full  interaction  kicks  and  whole  Keplerian  steps,  but  ending  with 
a  half  Kepler  step. 

This  ridiculously  simple  scheme  is  actually  a  remarkably  good  and  efficient 
integrator.  Wisdom  and  Holman  [10]  used  this  new  map  to  compute  the 
evolution  of  the  outer  planets  for  a  billion  years  and  compared  the  results  to 
the  results  of  the  845  million  year  integrations  performed  on  the  Digital  Orrery 
by  Sussman  and  Wisdom  [1]  using  conventional  integration  techniques.  All  of 
the  results  of  that  study  were  confirmed,  including  details  of  the  very  long- 
period  variations  (>  500  million  years)  in  the  inclination  of  Pluto,  and  the 
20  million  year  exponential  divergence  of  trajectories,  which  indicated  chaotic 
behavior  in  the  motion  of  Pluto.  The  new  map  is  about  an  order  of  magnitude 
faster  than  traditional  methods  of  integration. 

Our  100  Million  Year  Integration 

Our  100  million  year  integration  of  the  planetary  system  was  performed  us¬ 
ing  the  Supercomputer  Toolkit  [14].  The  Supercomputer  Toolkit  is  the  succes¬ 
sor  to  the  Digital  Orrery  [15].  It  is  a  small  multiprocessor  computer  optimized 
for  the  numerical  solution  of  systems  of  ordinary  differential  equations.  The 
Toolkit  was  built  as  a  collaboration  between  MIT  and  Hewlett-Packard.  It  is  a 
family  of  standard  software  and  hardware  modules  that  can  be  interconnected 
in  a  variety  of  configurations  as  appropriate  for  particular  applications.  Each 
processor  of  the  Toolkit  is  three  times  faster  than  the  entire  Distal  Orrery, 
as  measured  by  running  the  same  computation.  We  used  the  eight-processor 
configuration  of  the  Toolkit  at  MIT,  with  the  new  sjrmplectic  n-body  map,  to 
carry  out  eight  100  million  year  integrations  of  the  planetary  system.  Each 
processor  was  used  to  run  a  separate  integration  with  slightly  different  initial 
conditions  so  that  we  could  look  for  exponential  divergence. 

The  most  complete  of  the  long-term  integrations  of  the  whole  Solar  Sys¬ 
tem  is  that  of  QTD.  Particular  care  was  taken  by  QTD  to  make  their  physical 
model  accurate.  They  included  general  relativistic  corrections,  and  a  carefully 
crafted  quadrupole  approximation  to  account  for  the  long-term  effects  due  to 
the  finite  size  of  the  Earth-Moon  system.  They  used  initial  conditions  and 
masses  derived  from  JPL  ephemeris  DE102.  Our  physical  model  is  the  same 
as  that  of  QTD  except  in  our  treatment  of  the  effects  of  general  relativity. 
General  relativistic  corrections  can  be  written  in  Hamiltonian  form,  but  we 
have  not  been  able  to  integrate  them  analytically.  Instead  we  used  the  po¬ 
tential  approximation  of  Nobili  and  Roxburgh  [16],  which  is  easily  integrated, 
but  only  approximates  the  relativistic  corrections  to  the  secular  evolution  of 
the  shape  and  orientation  of  the  orbit. 

Multistep  methods  typically  require  more  than  a  hundred  steps  per  orbit  for 
stability;  the  mapping  method  is  stable  with  as  few  as  10  steps  per  orbit.  The 
step-size  for  our  integration  was  rather  arbitrarily  chosen  to  be  7.2  days;  with 
7.2-day  steps  output  points  are  easily  compared  to  the  ephemeris  of  QTD.  We 
integrated  backward  in  time.  A  22  million  year  integration  using  3.6-day  steps 
was  carried  out  to  check  that  the  7.2-day  integration  was  sufficiently  accurate. 
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All  of  the  position  calculations  were  done  in  pseudo-quadruple  precision,  in 
an  effort  to  reduce  the  effect  of  roundoff  error.  We  believe  now  that  the  use 
of  quadruple  precision  was  an  unnecessary  precaution,  and  needlessly  doubled 
the  computation  time  of  our  experiment.  Even  on  this  demanding  task  each 
processor  evolved  the  Solar  System  at  the  rate  of  about  30  years  per  second. 
Thus  our  eight  100  million  year  integrations  took  a  total  of  about  1000  hours 
of  Toolkit  time. 

We  have  compared  our  integration  to  the  reversed-time  segment  of  the 
0.5-day  step-size  integration  of  QTD,  and  the  agreement  is  quite  good.  The 
maximum  difference  in  the  argument  of  perihelion  of  Mercury  over  the  3  mil¬ 
lion  year  interval  is  of  order  0.0001  radians;  for  comparison,  the  precession  of 
the  argument  of  perihelion  due  to  general  relativity  is  about  2ir  radians  over 
3  million  years.  Evidently,  our  more  approximate  treatment  of  relativity  is  of 
little  consequence.  As  another  comparison.  Table  1  lists  the  maximum  differ¬ 
ences  between  the  eccentricities  of  the  planets  compared  to  QTD.  Also  listed 
sire  the  differences  between  the  integration  of  QTD  and  that  of  Laskar  [9].  The 


ILaskar  -  QTD|  [Toolkit  -  QTD| 


Mercury 

0.0041 

0.000018 

Venus 

0.0020 

0.000065 

Earth 

0.0024 

0.000059 

Mars 

0.0041 

0.000132 

Jupiter 

0.0038 

0.000047 

Saturn 

0.0081 

0.000162 

Uranus 

0.0051 

0.000008 

Neptune 

0.0026 

0.000002 

Pluto 

— 

0.000001 

Table  1:  The  maximum  differences  in  the  eccentricities  of  the  planets  in  the 
integrations  of  Laskar,  QTD,  and  the  Toolkit  show  excellent  agreement  among 
the  integrations.  Laskar  did  not  include  Pluto  in  his  calculation. 

table  illustrates  that  the  evolution  computed  vrith  the  mapping  agrees  qidte 
well  with  the  more  conventional  direct  integration  of  QTD.  The  Toolkit  inte¬ 
gration  and  QTD  are  mutually  more  consistent  than  either  is  to  the  integration 
of  Laskar,  though  it  is  not  clear  whether  this  discrepancy  is  due  primarily  to 
model  differences  or  to  the  approximations  used  by  Laskar.  The  model  dif¬ 
ferences  are  actually  to  our  benefit  because  they  help  address  the  important 
question  of  the  sensitivity  of  our  results  to  model  parameters. 

The  Evolution  of  the  Solar  System  is  Chaotic 

Exponential  divergence  of  nearby  trajectories  is  indicative  of  chaotic  behav¬ 
ior.  The  divergence  of  trajectories  may  have  both  quasiperiodic  and  exponen¬ 
tial  components.  Of  course,  if  present,  the  exponential  component  eventually 
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dominates.  In  order  to  estimate  the  Lyapunov  exponent,  we  extracted  the  ex¬ 
ponential  part  of  the  divergence  from  the  quasiperiodic  part  by  the  following 
procedure.  We  first  converted  the  positions  and  momenta  into  Keplerian  ele¬ 
ments,  then  formed  the  standard  shape  and  orientation  variables  h  =  e  sin  to, 
k  =  e  cos  tu,  p  =  sin  i/2  sin  H,  and  q  =  sin  t/2  cos  fl,  where  e  is  the  eccentricity, 
w  is  the  longitude  of  perihelion,  t  is  the  inclination,  and  SI  is  the  longitude 
of  the  ascending  node.  The  full  set  of  these  variables  for  all  of  the  planets 
constitutes  the  *Tull  secular  phase  space.”  Distance  in  this  space  is  the  or¬ 
dinary  Euclidean  distance.  The  quasiperiodic  component  is  not  as  strong  in 
the  secular  phase  space  as  it  is  in  the  full  phase  space,  allowing  us  to  study 
the  exponential  part  of  the  divergence  more  easily.  If  exponential  divergence  is 
observed  in  the  secular  phase  space  exponential  divergence  will  also  be  present 
in  the  full  phase  space. 

The  divergence  in  the  secular  phase  space  between  two  of  the  100  million 
year  calculations  is  shown  in  Figure  1.  The  initial  conditions  for  the  two  calcu¬ 
lations  differed  by  about  1  millimeter  in  the  x  coordinate  of  Pluto.  The  secular 
divergence  has  two  distinct  segments.  In  the  latter  segment  the  divergence  is 
dominated  by  an  exponential  with  a  timescale  of  about  4  million  years.  The 
initial  segment  is  apparently  dominated  by  a  smaller  exponent  with  a  timescale 
of  about  12  million  years. 

If  the  system  is  chaotic  then  the  divergence  of  individual  planets  will  also 
be  exponential  with  the  same  exponent  as  the  whole  system  if  the  trajectories 
are  followed  for  sufficiently  long  time.  Over  shorter  intervals  the  divergence 
of  individual  planets  can  be  different,  and  examination  of  the  individual  di¬ 
vergences  can  give  insight  into  the  mechanisms  responsible  for  the  chaotic 
behavior  in  the  system. 

In  the  inner  Solar  System  the  individual  planet  divergences  in  the  secu¬ 
lar  phase  space  are  similar  to  that  of  the  full  secular  phase  space,  with  two 
distinct  segments.  In  the  outer  Solar  System,  over  most  of  the  100  million 
years  spanned  by  our  integrations  the  divergences  appear  to  be  dominated 
by  a  12  million  year  exponential  component,  with  evidence  of  the  4  million 
year  component  appearing  in  only  the  last  5  million  years.  A  change  in  the 
timescale  of  exponential  divergence  could  occiir  if  the  system  went  from  one 
region  of  the  phase  space  characterized  by  one  exponent  into  another  region 
characterized  by  the  other  exponent.  This  possibility  can  be  ruled  out  because 
the  4  million  year  component  appears  much  later  in  the  outer  planets  than  in 
the  inner  planets.  A  viable  interpretation  is  that  there  are  two  distinct  mech¬ 
anisms  generating  exponential  divergence  that  are  simultaneously  operating. 
Apparently  the  inner  planets  are  more  sensitive  indicators  of  the  4  million  year 
process  than  are  the  outer  planets. 

The  evolution  of  Pluto  over  the  full  100  million  years  has  characteristics 
quite  similar  to  those  found  in  long  term  integrations  of  the  outer  planets. 
For  example,  we  observe  the  34  million  year  amplitude  modulation  of  the  3.8 
million  year  oscillation  of  the  argument  of  perihelion  of  Pluto  [5]. 

Figure  2  shows  the  exponential  divergence  of  the  difference  of  positions  of 


5 


Figure  1;  Exponential  divergence  of  nearby  trajectories  is  indicated  by  the 
average  linear  growth  of  the  logarithm  of  the  full  secular  phase  space  distance 
between  two  different  runs  with  very  slight  differences  in  initial  conditions.  The 
segment  following  the  initial  transient  has  an  exponential  timescale  of  about 
12  million  years.  The  divergence  is  subsequently  dominated  by  an  exponential 
with  a  timescale  near  4  million  years. 
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time  [10*  years] 


Figure  2:  The  diveq^ce  of  the  distance  between  the  positions  of  Pinto  in  two 
different  runs  with  very  slight  differences  in  initial  conditions  is  characterized 
by  a  timescale  near  12  million  years. 


Pluto  in  two  of  our  100  million  year  runs.  These  runs  differed  by  1  part  in  10^® 
in  one  coordinate  of  the  initial  position  of  Mars.  The  plot  is  consistent  with 
an  exponential  divergence  of  about  12  million  years.  This  rate  of  exponential 
divergence  is  similar  to  that  observed  in  the  outer  planets,  and  to  the  slower 
component  observed  in  the  full  secular  phase  space.  One  interpretation  is 
that  Pluto  is  passively  driven  by  chaotic  behavior  of  the  rest  of  the  Solar 
System.  Another  interpretation  is  that  Pluto  is  independently  chaotic  with  a 
rate  of  exponential  divergence  which  is  only  coincidentally  similar  to  that  of 
the  rest  of  the  system.  We  present  evidence  for  the  latter  interpretation  below. 
Considering  the  typical  slow  convergence  of  Lyapunov  estimates,  a  12  million 
year  timescale  of  exponential  divergence  is  in  satisfactory  agreement  with  the 
inverse  Lyapunov  exponent  of  15  to  20  million  years  for  Pluto  reported  by 
Sussman  and  Wisdom  [1]. 

Secular  Resonances 

Laskar  [17,  18]  has  found  three  resonance  angles  which  alternately  librate 
and  circulate.  This  may  be  an  important  corroboration  of  the  chaos  as  indi¬ 
cated  by  the  exponential  divergence.  The  Laskar  angles  are 

=  K  -  o  -  («;  -  ftS).  (3) 

=  2iwl  -  wl)  -  (ftj  -  ill),  (4) 

^3  =  (mj  -  -  (ftj  -  ft5).  (5) 

The  quantities  and  flf  are  the  angles  of  the  proper  modes  as  defined  by 
Laskar  [17].  On  the  basis  of  the  alternate  libration  of  the  angles  o’}  and  <^3, 
which  cannot  simultaneously  librate,  Laskar  concluded  that  the  mechanism 
responsible  for  the  chaotic  behavior  of  the  Solar  System  was  resonance  overlap 
of  the  corresponding  secular  resonances. 

In  our  calculation,  <ri  and  also  alternately  circulate  and  librate,  but  <73 
just  circulates,  as  shoum  in  Figure  3.  Our  angles  track  Laskar’s  angles  for  the 
initial  segment  of  the  computations,  but  they  soon  diverge  enough  so  that  the 
changes  between  libration  and  drctilation  occur  at  very  different  times.  Such 
differences  are  consistent  with  the  chaotic  character  of  the  evolution.  The 
higher-order  angle 

<r4  =  3K-w|)-2(nS-ft5)  (6) 

is  also  incompatible  with  <73  and  <7,.  Xn  our  calculation  <74  has  intervals  of 
libration.  The  data  presented  by  Laskar  suggest  that  this  is  also  the  case  in 
his  calculation,  though  the  angle  is  not  explicitly  mentioned.  Since  1/1  < 
3/2  <  2/1  it  is  possible  that  we  and  Laskar  are  seeing  different  portions  of 
a  chaotic  zone  that  spans  the  phase  space  from  the  1/1  resonance  to  the  2/1 
resonance.  In  addition  to  these  angles  we  also  find  that  the  angle 

<73  =  (m?  -  mS)  -»-  (ft;  -  ftj)  (7) 


alternates  between  circulation  and  libration. 


8 


time  [10®  years] 


Figure  3:  In  our  calculation  two  of  the  angles  of  Laskar,  <ri  and  o’],  alternately 
circulate  and  librate.  Laskar’s  third  angle,  o’],  circulates.  The  higher  order 
angle  o’4,  which  is  incompatible  with  both  o’]  and  O’],  has  intervals  of  libration. 
The  angle  o']  also  alternately  circulates  and  librates. 


Though  our  calculations  c^e  completely  consistent  with  those  of  Laskar,  we 
are  not  fully  convinced  that  his  proposed  mechanism  accounts  for  the  observed 
chaotic  behavior  of  the  Solar  System.  First,  it  is  not  clear  that  the  alternate 
libration  and  circulation  of  is  indicative  of  a  dynamically  significant  chaotic 
separatrix.  In  Figure  4  we  show  a  polar  plot  of  the  amplitude  and  phase  of 
the  combination  of  proper  modes  corresponding  to  ai  and  <73.  Our  polar  plot 
of  (Ti  is  similar  to  that  shown  in  Laskar  [17].  However,  our  polar  plot  of  <73 
differs  from  that  shown  in  Laskar  [17]  and  Laskar,  Quinn,  and  Tremaine  [9]. 
For  the  polar  amplitude  they  use  just  the  amplitude  of  eccentricity  mode  4, 
which  is  never  close  to  zero.  The  amplitude  corresponding  to  <73  has  additional 
factors,  some  of  which  approach  zero.  Using  just  the  amplitude  of  eccentricity 
mode  4  gives  an  artificial  impression  of  a  transition  from  libration  to  -ircula- 
tion.  Our  plot,  which  includes  all  the  factors  in  the  amplitude,  does  not  give 
a  clear  indication  of  a  chaotic  separatrix.  Rather,  it  gives  the  impression  of 
a  complex  high-dimensional  trajectory  projected  onto  a  plane  (see  Figure  4). 
The  alternation  of  libration  and  circulation  of  the  polar  angle  may  just  be 
an  artifact  of  the  origin  being  in  the  midst  of  this  complex  projected  trajec¬ 
tory.  A  simple  integrable  model  for  motion  near  the  2/1  commensurability  in 
the  planar-elliptic  restricted  three-body  problem  [19,  20]  illustrates  how  the 
alternate  drctdation  and  libration  of  an  angle  can  be  misleading.  There  the 
resonance  angle  corresponding  to  the  largest  term  in  the  disturbing  function 
can  show  alternate  circulation  and  libration,  even  though  there  is  no  chaotic 
behavior.  The  phenomenon  can  be  eliminated  by  an  appropriate  change  of 
variables.  The  polar  plot  corresponding  to  <ri  is  more  like  that  expected  of 
an  angle  associated  with  a  chaotic  mechanism.  It  not  only  alternately  circu¬ 
lates  and  librates,  but  as  it  circulates  it  loiters  around  an  apparent  unstable 
equilibrium.  There  is  also  an  excluded  re^on  near  the  center  of  the  plot. 

Furthermore,  there  are  too  many  unrelated  angles  which  alternately  cir- 
ctilate  and  librate.  It  might  be  expected  that  only  one  incompatible  set  of 
angles  would  show  alternate  circulation  and  libration.  However,  we  see  the 
phenomenon  in  angles  involving  unrelated  modes.  Eccentricity  and  inclina¬ 
tion  modes  3  and  4  are  involved  in  <73,  <73,  and  474;  eccentricity  modes  1,  5,  and 
8,  and  inclination  modes  1,  2,  and  8  are  involved  in  <ri  and  475.  The  two  sets 
of  modes  are  disjoint,  and  yet  there  are  correlations  in  the  behavior  of  angles 
associated  with  unrelated  sets  of  modes.  Most  striking  is  the  transition  in  be¬ 
havior  in  four  of  the  angles  near  an  integration  time  of  50  million  years.  This 
suggests  that  a  single  mechanism  is  driving  all  of  the  angles.  If  the  mechanism 
is  associated  with  one  of  the  angles  presented,  we  feel  the  most  convincing  is 
<7}.  It  is  also  possible  that  the  mechanism  generating  the  chaos  is  unrelated 
to  all  of  these  angles,  and  that  they  are  all  just  sensitive  indicators  of  chaotic 
irregularity  of  the  underlying  system  trajectory. 

We  have  not  Mentified  any  other  angles  whose  motion  is  suggestive  of  a 
dynamical  mechanism.  The  LONGSTOP  team  speculated  [21]  that  the  angle 
2wl  —  2c77  +  07  —  fig  might  alternately  circulate  and  librate,  and  thereby 
provide  evidence  of  chaotic  behavior  of  the  outer  planets.  We  find  that  this 
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Figure  4:  The  variables  are  Xj  =  Ai  cosffi  and  y,-  =  sin  where  Ai  are  the 
amplitudes  of  the  combination  of  the  Laskar  proper  modes  corresponding  to 
the  angles  The  plot  corresponding  to  <73  gives  the  impression  of  a  complex 
high-dimensional  trajectory  projected  onto  a  plane.  The  plot  corresponding 
to  0*1  shows  some  indication  of  a  dynamically  significant  chaotic  sepaiatrix. 
There  is  a  suggestion  of  an  unstable  equilibrium  and  associated  asymptotic 
trajectories,  and  there  is  a  region  near  the  center  that  is  avoided. 


angle  rotates  uniformly  with  a  period  somewhat  greater  than  50  million  years. 
This  angle  is  one  of  30  angles  presented  by  the  LONGSTOP  team  in  their 
“numerology  table.”  We  find  all  30  angles  rotate  uniformly;  there  is  no  hint 
of  chaotic  irregularity  in  the  behavior  of  any  of  the  angles. 

That  particular  secular  resonances  are  responsible  for  the  chaotic  behavior 
of  the  Solar  System  could  be  established  by  an  analytical  demonstration  that 
the  resonances  involved  are  sufficiently  strong  and  close  to  give  resonance 
overlap.  At  the  moment,  however,  we  feel  no  dynamical  mechanism  for  the 
observed  chaotic  behavior  of  the  Solar  System  has  been  clearly  demonstrated. 

Chaotic  Evolution  of  the  Jovian  Planets 

Despite  the  growing  number  of  long-term  integrations  of  the  outer  plan¬ 
ets,  to  date  no  integration  has  directly  tested  whether  the  massive  planets 
themselves  evolve  chaotically.  In  our  845  million  year  integration  of  the  outer 
planets  [1],  the  orbital  elements  of  Neptune  appeared  to  have  discrete  line  spec¬ 
tra,  in  marked  contrast  to  the  clearly  broad-band  character  of  the  spectrum 
of  Pluto.  On  the  basis  of  this  spectral  evidence  we  dismissed  the  possibility 
of  chaotic  behavior  in  the  Jovian  planets.  To  be  thorough,  we  have  now  car¬ 
ried  out  a  new  billion  year  evolution  of  the  outer  planets,  using  the  mapping 
method,  with  a  slightly  perturbed  initial  position  of  Neptune.  We  found,  to  our 
surprise,  that  the  subsequent  evolution  of  the  Jovian  planets  diverged  expo¬ 
nentially  from  the  first  calculation,  with  a  timescale  of  exponential  divergence 
of  only  5  million  years! 

Our  initial  reaction  was  that  there  must  be  something  wrong  with  our 
method  of  integration.  To  check  this  we  carried  out  two  new  integrations 
of  the  outer  planets  spanning  100  million  years  using  the  traditional  linear 
multistep  Stormer  predictor,  the  same  integrator  used  in  the  Digital  Orrery 
integrations.  The  initial  conditions  and  masses  were  the  same  as  in  the  Dig¬ 
ital  Orrery  integrations.  The  integrations  were  carried  out  in  ordinztry  IEEE 
double  precision  (64  bits)  with  a  step-size  of  32.7  days,  the  same  step-size 
used  in  the  845  million  year  Digital  Orrery  integrations.  We  found  that  the 
trajectories  of  the  Jovian  planets  diverged  exponentially  with  a  timescale  of 
about  19  million  years.  Apparently,  we  were  misled  by  the  spectral  evidence. 

To  further  check  that  this  result  did  not  depend  on  either  the  step-size 
or  the  precision  of  the  calculation  we  carried  out  four  more  integrations  of 
the  outer  planets  using  the  Stormer  predictor.  Each  spanned  more  than  400 
million  years.  In  these  integrations  the  accumulation  of  position  was  carried 
out  in  pseudo-quadruple  precision,  as  in  the  Digital  Orrery  integrations.  One 
pair  used  a  step-size  of  32.7  days,  the  special  Digital  Orrery  step-size;  the  other 
pair  used  an  arbitrarily  chosen  step-size  of  28  days.  The  initial  conditions  were 
the  same  as  in  our  earlier  outer  planet  integrations.  In  one  run  of  each  pair 
the  initial  position  of  Neptune  was  perturbed  by  7.5mm.  The  energy  errors 
again  grew  linearly  with  time,  with  slopes  between  that  of  our  210  million  year 
Digital  Orrery  integration  and  that  of  our  845  million  year  integration.  Both 
pairs  of  runs  gave  remarkably  consistent  restilts.  The  secular  divergence  of  the 
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Jupiters  in  the  28-day  runs  is  shown  in  Figure  5.  The  Jovian  planets  diverged 
exponentially  with  a  timescale  somewhat  longer  than  20  million  years. 

The  map  and  Stormer  calculations  both  indicate  that  the  motion  of  the 
Jovian  planets  is  chaotic.  However,  they  are  discrepant  in  the  estimate  of 
the  exponential  timescale.  We  believe  this  is  a  resiilt  of  the  fact  that  the 
trajectory  computed  by  the  mapping  method  is  not  exactly  the  same  as  the 
trajectory  determined  by  the  initial  conditions.  The  mapping  differs  from  the 
actual  n-body  dynamics  by  the  addition  of  extra  high-frequency  terms.  The 
averaging  principle  suggests  that  these  high-frequency  terms  do  not  contribute 
to  the  long-term  evolution,  and  the  close  agreement  of  our  results  with  those  of 
Laskar  once  again  confirms  the  validity  of  the  use  of  averaging.  However,  the 
initial  conditions  used  in  an  averaged  system  should  properly  take  into  account 
the  presence  of  the  extra  high-frequency  components  in  the  unaveraged  system. 
The  use  of  the  same  initial  conditions  in  the  averaged  and  unaveraged  systems 
corresponds  to  slightly  different  initial  conditions  for  the  long-term  evolutions. 
With  the  mapping  we  do  not  yet  know  how  to  properly  adjust  the  initial 
conditions  to  account  for  the  extra  high-frequency  terms.  In  our  outer-planet 
integrations  these  slight  uncontrolled  adjustments  of  effective  initial  conditions 
appear  to  yield  different  estimates  of  the  Lyapunov  exponent. 

To  investigate  this  further  we  carried  out  8  integrations  of  the  outer  plan¬ 
ets  using  the  map  with  different  step-sizes,  ranging  from  0.3  years  to  1  year. 
Changing  the  step-size  changes  the  high-frequency  components,  and  slightly 
changes  the  effective  initial  condition  for  the  long-term  evolution.  Each  inte¬ 
gration  spanned  about  300  million  years.  We  found  that  the  measured  diver¬ 
gence  timescale  varied  from  about  3  million  years  to  about  30  million  years. 
The  dispersion  in  the  estimates  of  the  Lyapunov  exponent  are  much  larger 
than  the  dispersion  observed  in  the  Stormer  runs.  The  Lyapunov  exponent 
was  not  obviously  correlated  with  step-size,  in  particular  the  estimate  of  the 
Lyapunov  exponent  was  not  monotonic  with  step-size.  In  one  of  these  runs, 
with  step-size  near  0.617979  years,  the  motion  of  the  outer  planets  was  clearly 
quasiperiodic;  the  secular  phase  space  divergence  did  not  grow  exponentially. 
The  results  suggest  that  the  Lyapunov  exponent  for  the  Jovian  planets  is  not 
a  simple  function  of  the  initial  conditions.  Most  nearby  initial  conditions  lead 
to  exponential  divergence  (most  with  a  shorter  timescale  for  exponential  di¬ 
vergence  than  that  obtained  with  the  Stormer  integrations),  but  there  are  also 
nearby  initial  conditions  that  do  not  give  chaotic  behavior. 

With  each  outer  planet  integration  we  ran  a  pair  of  massless  Plutos,  with 
initial  positions  differing  by  about  1cm.  A  remarkable  result  is  that  the  expo¬ 
nential  divergence  of  Plutos  always  has  a  timescale  between  10  and  20  million 
years,  independent  of  how  chaotically  the  Jovian  planets  behave.  This  is  true 
even  in  the  most  extreme  runs,  where  the  Jovian  planets  were  quasiperiodic, 
and  where  the  Jovian  planets  diverged  with  a  timescale  of  3  million  years.  This 
clearly  demonstrates  that  the  mechanism  generating  the  chaotic  behavior  in 
the  motion  of  Pluto  is  extremely  robust,  and  independent  of  whether  the  rest 
of  the  system  is  chaotic. 
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Figure  5:  The  secular  phase  space  divergence  of  the  trajectories  of  Jupiter  in 
the  28  day  Stormer  integrations  show  a  timescale  for  exponential  divergence 
that  is  somewhat  longer  than  20  million  years.  The  divergence  saturates  after 
about  250  million  years  at  a  rather  small  value,  perhaps  indicating  a  hidden 
constraint  on  the  trajectories. 


Numerical  Chaos? 

The  fact  that  almost  all  long-term  integrations  of  the  Solar  System  give 
exponential  divergence  of  trajectories  with  a  timescale  in  the  range  of  3-30 
million  years  in  physically  quite  different  models  is  very  striking  and  unsettling. 
Another  unsettling  feature  of  the  chaotic  behavior  we  observe  in  long-term 
planetary  integrations  is  that  nothing  dramatic  happens.  This  is  compounded 
by  the  fact  that  in  no  case  have  we  unambiguously  identified  the  mechanism 
responsible  for  the  chaotic  behavior.  This  lack  of  mechanism  and  lack  of 
obvious  irregular  behavior  is  in  marked  contrast  to  the  clearly  understood 
mechanisms  and  irregular  character  observed  in  other  examples  of  chaotic 
behavior  in  the  solar  system,  for  example,  asteroids  on  chaotic  trajectories 
near  commensurabilities  [11,  12,  22],  the  chaotic  tumbling  of  Hyperion  [23] 
and  other  irregularly  shaped  satellites  [24],  and  the  chaotic  motion  of  the 
Uranian  satellites  near  commensurabilities  [25,  26,  27]. 

Perhaps  the  exponential  divergence  is  a  numerical  artifact?  The  detailed 
agreement  of  the  billion  year  evolution  of  Pluto  using  different  integrators  is 
impressive  evidence  against  this  perverse  possibility.  Furthermore,  the  detailed 
agreement  between  our  100  million  year  Solar  System  integration  and  that  of 
Laskar  is  particularly  convincing,  because  of  the  radically  different  methods 
used. 

To  further  convince  ourselves  that  not  all  long-term  integrations  are  subject 
to  some  universal  numerical  instability,  which  yields  a  spurious  exponential  di¬ 
vergence,  we  have  carried  out  a  control  integration.  We  have  integrated  the 
outer  planets  without  Uranus  for  about  250  million  years,  with  the  mapping 
method.  Over  that  period  we  see  no  sign  of  exponential  divergence  of  the  re¬ 
maining  Jovian  planets.  This  integration,  together  with  the  isolated  quasiperi- 
odic  integration  mentioned  above,  shows  that  numerical  models  of  planetary 
systems  can,  in  fact,  show  quasiperiodic  behavior  on  a  several  hundred  million 
year  time-scale.  Long-term  integrations  do  not  always  give  positive  Lyapunov 
exponents. 

Altogether,  the  evidence  for  the  chaotic  behavior  in  these  long-term  plan¬ 
etary  integrations  is  very  convincing,  but  there  remains  the  logical  possibility 
that  the  exponential  divergence  is  a  subtle  numerical  artifact.  To  positively 
conclude  that  the  chaos  observed  in  these  long-term  planetary  integrations 
is  not  a  result  of  numerical  artifacts  requires  an  unambiguous  identification 
of  a  physical  mechanism  and  an  analytic  evaluation  to  determine  that  the 
mechanism  actually  accounts  for  the  observed  chaos. 

Conclusions  and  Speculations 

Our  100  million  year  integration  of  the  entire  Solar  System  indicates  that 
the  Solar  System  is  chaotic  with  a  timescale  for  exponential  divergence  of 
about  4  million  years.  The  fact  that  we  find  similar  behavior  in  all  respects 
to  the  calculation  of  Laskar  strongly  supports  the  conclusion  that  the  Solar 
System  is  chaotic.  That  we  and  Laskar  have  carried  out  different  kinds  of 
numerical  experiments,  with  slightly  different  masses,  slightly  different  initial 
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conditions,  and  slightly  different  physics,  shows  that  the  chaotic  character  of 
the  Solar  System  is  not  sensitively  dependent  on  the  precise  model  or  numerical 
methods. 

Our  experiments  indicate  that  the  Jovian  planets  by  themselves  behave 
chaotically  for  most  initial  conditions  near  our  reference  system,  though  our 
estimates  of  the  Lyapunov  exponent  have  rather  large  dispersion. 

All  of  our  estimates  of  the  Lyapunov  exponent  of  Pluto  ^ve  approximately 
the  same  divergence  timescale  of  10*20  million  years,  with  different  methods 
of  integration,  different  planetary  masses,  different  initial  conditions,  and  even 
independently  of  whether  the  rest  of  the  system  is  behaving  chaotically.  Our 
earlier  result  that  the  evolution  of  Pluto  is  chaotic  is  thus  multiply  confirmed. 

We  will  not  fully  understand  the  consequences  of  the  observed  chaotic  evo¬ 
lution  of  the  Solar  System  until  we  clearly  understand  the  dynamical  mecha¬ 
nisms  responsible  for  it.  Though  in  our  calculation  the  behavior  of  the  secular 
resonance  angles  is  consistent  with  those  in  Laslcar’s  calculation,  the  identifi¬ 
cation  of  resonance  overlap  of  particular  secular  resonances  as  the  mechanism 
generating  the  chaotic  behavior  of  the  Solar  System  is  not  unambiguously 
demonstrated.  Our  numerical  experiments  suggest  that  there  are  at  least  two 
independent  mechanisms  generating  chaotic  behavior.  One  mechanism  gen¬ 
erates  chaos  in  the  Jovian  subsystem,  and  Pluto  is  independently  chaotic  in 
the  field  of  the  Jovian  planets.  Yet  another  mechanism  is  suggested  by  the 
simultaneous  presence  of  two  different  exponential  timescales  in  our  full  So¬ 
lar  System  integrations.  Secular  resonances  among  the  inner  planets  may  be 
driving  the  faster  timescale,  as  Laskar  suggested.  However,  the  most  convinc¬ 
ing  of  the  secular  resonances  involve  Mercury.  Both  Mercury  and  Pluto  have 
high  eccentricity  and  inclination,  which  strongly  couples  the  eccentricity  and 
inclination  secular  subsystems.  Perhaps  one  of  the  mechanisms  generating  the 
chaos  originates  with  Mercury,  and  is  similar  to  the  mechanism  generating  the 
chaos  in  the  motion  of  Pluto? 

References 

[1]  G.J.  Sussman  and  J.  Wisdom  Science  241,  433  (1988). 

[2]  J.  Laskar  Nature  3S8,  237  (1989). 

[3]  C.J.  Cohen,  E.C.  Hubbard,  and  C.  Oestermnter,  in  Astronomical  Papers 
of  the  American  Epbemeris  (Government  Printing  Office,  Washington, 
DC,  1973),  vol.  22,  pp.  1-92. 

[4]  H.  Kinoshita,  and  H.  Nakai  Celestial  Mechanics  S4,  203  (1984). 

[5]  J.  Applegate,  M.R.  Douglas,  Y.  Gursel,  G.J.  Sussman,  and  J.  Wisdom 
Astron.  J.,  02,  176  (1986). 


16 


[6]  A.E.  Roy,  I.W.  Walker,  A.J.  Macdonald,  I.P.  WiUiams,  K.  Fox,  C.D. 
Murray,  A.  Milani,  A.M.  Nobili,  P.J.  Message,  A.T.  Sinclair,  and  M. 
Carpino,  Vistas  in  Astronomy  S2,  95  (1988). 

[7]  D.L.  Richardson,  and  G.F.  Walker,  in  Astrodynamics  1987  (Advances 
in  the  Astronautical  Sciences,  Vol.  65),  J.K.  Soldner,  A.K.  h^ra,  R.E. 
Lindberg,  and  W.  Williamson,  eds.  (Univelt,  San  Diego,  1987),  1473. 

[8]  T.R.  Quinn,  S.D.  Tremaine,  and  M.  Duncan  Astron.  3. 101,  2287  (1991). 

[9]  J.  Laskar,  T.R.  Quinn,  and  S.D.  Tremaine,  Icarus  05,  148  (1992). 

[10]  J.  Wisdom,  and  M.  Holman  Astron.  J.  102, 1528  (1991). 

[11]  J.  Wisdom  Astron.  J.  87,  577  (1982). 

[12]  J.  Wisdom  Icarus  56,  51  (1983). 

[13]  V.I.  Arnold  Mathematical  Methods  of  Classical  Mechanics  (Springer- 
Verlag,  1978). 

[14]  H.  Abelson,  A.A.  Berlin,  J.  Katzenelson,  W.  McAllister,  G.J.  Rozas,  G.J. 
Sussman,  and  J.  Wisdom,  preprint. 

[15]  J.  Applegate,  M.R.  Douglas,  Y.  Gursel,  P.  Hunter,  C.  Seitz,  and  G.J. 
Sussman,  IEEE  Trans.  Comput.  CS4(1985). 

[16]  A.  Nobili  and  I.  Roxburgh,  in  Relativity  in  Celestial  Mechanics  and  As¬ 
trometry,  Kovalevsky,  Brumberg  eds.  (Reidel,  Dordrecht,  1986). 

[17]  J.  Laskar  Icarus  88,  266  (1990). 

[18]  J.  Laskar  in  Proc.  lAU  Symposium  152^  Ferraz-Mello  ed.  (1991). 

[19]  J.  Wisdom  Celestial  Mechanics  58,  175  (1986). 

[20]  J.  Henrard,  A.  Lemaitre,  A.  Milani,  and  C.D.  Murray  Celestial  Mechanics 
58,  335  (1986). 

[21]  A.  Nobili,  A.  Milani,  and  M.  Carpino  Astron.  Astrophys.  210,  313  (1989). 

[22]  J.  Wisdom  Icarus,  65,  272  (1985). 

[23]  J.  Wisdom,  S.J.  Peale,  and  F.  Mignard  Icarus  58, 137  (1984). 

[24]  J.  Wisdom  Astron.  J.  04,  1350. 

[25]  W.  Tittemore  and  J.  Wisdom  Icarus  74,  172  (1988). 

[26]  W.  Tittemore  and  J.  Wisdom  Icarus  78,  63  (1989). 

[27]  W.  Tittemore  and  J.  Wisdom  Icarus  85,  394  (1990). 


17 


[28]  This  work  depended  upon  many  people.  The  Supercomputer  Toolkit  com¬ 
puter  was  designed  and  built  with  the  help  and  generous  support  of  the 
Hewlett  Packard  Company.  We  especially  thank  Joel  Bimbaum  of  Hewlett 
Packard  for  his  personal  support.  William  McAllister  led  the  team  at  HP. 
Dan  Zuras  of  HP  helped  with  the  development  of  the  scientific  subrou¬ 
tines  we  used.  At  MIT  Andy  Berlin,  Bill  Rozas,  Henry  Wu,  Hal  Abelson, 
and  Jacob  Katzenelson  contributed  to  the  hardware  and  software  de¬ 
sign.  Matthew  Holman  helped  with  the  development  and  testing  of  the 
mapping  integrator.  We  thank  Scott  Tremaine,  Jacques  Laskar,  Thomas 
Quinn,  and  Martin  Duncan  for  hdpful  discussions. 

This  report  describes  research  done  at  the  Artificial  Intelligence  Labora¬ 
tory  of  the  Massachusetts  Institute  of  Technology.  Support  for  the  labora¬ 
tory’s  artificial  intelligence  research  is  provided  in  part  by  the  Advanced 
Research  Projects  Agency  of  the  Department  of  Defense  under  Office  of 
Naval  Research  contract  NOOO 14-89- J-3202  and  in  part  by  the  National 
Science  Foundation  under  grant  number  MIP-9001651.  Jack  Wisdom  is 
also  supported  in  part  by  a  NSF  Presidential  Young  Investigator  Award 
AST-887365  and  in  part  by  the  NASA  Planetary  Geology  and  Geophysics 
Program  under  grant  NAGW-706.  We  are  grateful  for  the  hospitality  of 
the  California  Institute  of  Technology  where  both  authors  are  on  leave 
from  MIT. 


18 


